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Abstract 

In this paper the space of images is considered as a Riemannian manifold using the metamorphosis ap¬ 
proach IM Y01 1 It Y 05al ITY05bl . where the underlying Riemannian metric simultaneously measures the cost 
of image transport and intensity variation. A robust and effective variational time discretization of geodesics 
paths is proposed. This requires to minimize a discrete path energy consisting of a sum of consecutive image 
matching functionals over a set of image intensity maps and pairwise matching deformations. For square- 
integrable input images the existence of discrete, connecting geodesic paths defined as minimizers of this 
variational problem is shown. Furthermore, F-convergence of the underlying discrete path energy to the con¬ 
tinuous path energy is proved. This includes a diffeomorphism property for the induced transport and the 
existence of a square-integrable weak material derivative in space and time. A spatial discretization via fi¬ 
nite elements combined with an alternating descent scheme in the set of image intensity maps and the set 
of matching deformations is presented to approximate discrete geodesic paths numerically. Computational 
results underline the efficiency of the proposed approach and demonstrate important qualitative properties. 


1 Introduction 

The study of spaces of shapes from the perspective of a Riemannian manifold allows to transfer many important 
concepts from classical geometry to these usually infinite-dimensional spaces. During the past decade, this 
Riemannian approach had an increasing impact on the development of new methods in computer vision and 
imaging, ranging from shape morphing and modeling, e.g. OKMP07I . and shape statistics, e.g. IIFLPJ04I . to 
computational anatomy IIBMTY02I . A variety of Riemannian shape spaces has been investigated in the liter¬ 
ature. Some of them are finite-dimensional and consider polygonal curves or triangulated surfaces as shapes 
IIKMP071ILSDM101 . but most approaches deal with infinite-dimensional spaces of shapes. Prominent examples 
with a full-fledged geometric theory are spaces of planar curves with curvature-based metric IIMM06I . elastic 
metric IISJJK06I or Sobolev-type metric IICKPF05IIMM07IISYM07I . The concept of optimal transport was 
used to study the space of images, where image intensity functions are considered as probability measures, 
e.g. Zhang et al. IIZYHT07I minimize the Monge-Kantorovich functional Jjj\iIj(x) — x\'^ po{x) dx over all mass 
preserving mappings -0 Benamou and Brenier BBBOOI used a flow reformulation of optimal transport, 

which nicely fits into the Riemannian context. 

For only a few nontrivial application-oriented Riemannian spaces geodesic paths can be computed in closed 
form (e.g. IYMSM08I ISMSYl fl '). else the system of geodesic ODEs has to be solved using numerical time 
stepping schemes (e.g. IIKSMJ04I IBMTYOSl ). Alternatively, geodesic paths connecting shapes can also be 
approximated via the minimization of discretized path length IISCC06II or path energy functionals IIFJSY09I 
IWBRSl fl . In this paper, we will develop such a variational time discretization on the space of images using 
the metamorphosis approach proposed by Trouve and Younes BTYOSbl lTY05al [HTY09I . This approach is a 
generalization of the flow of diffeomorphism approach initiated by Dupuis, Grenander and Miller IIDGM98I . 

The concept of variational time discretization is a powerful tool in the discretization of gradient flows and 
for Hamiltonian mechanical systems. The analog of the time discrete path energy considered here is a dis¬ 
crete action sum. For a historic account we refer to IIHLW06I . Numerical analysis was exploited from the 
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F-convergence perspective in IIMO04I . and from the ODE-discretization perspective under the name of varia¬ 
tional integrators in BLMOW041 lOB JM 111 . Thereby, the time continuous Lagrangian on some time interval is 
replaced by a time discrete functional related to our functional W and defined directly on configuration variables 
and not involving momentum variables. 

Instead of discretizing the underlying flow and incorporating the target configuration at the end time via a 
constraint, the variational discretization is based on the direct minimization of a discrete path energy subject to 
data given at the initial and the end time. This approach turned out to be very stable and robust, and even for 
very small numbers of time steps one obtains qualitatively good results. Furthermore, proceeding from coarse 
to fine time discretization, an efficient cascadic minimization strategy can be implemented. In the context of 
shape spaces, this concept has already been used in the space of viscous objects OFJSY09I IWBRSlTl IRW131 . 
but without a rigorous mathematical foundation. In IIRW14I . a discrete geodesic calculus on finite- and on 
certain infinite-dimensional shape spaces with the structure of a Hilbert manifolds was developed and a full- 
fledged convergence analysis could be established. This theory immediately applies for instance to the (finite- 
dimensional) Riemannian manifold of discrete shells llHRWW12llHRS~*~14i . In this paper, we expand part of 
this theory to the metamorphosis model, which lacks a Hilbert manifold structure. 

In what follows, we will briefly review both the flow of diffeomorphism and the metamorphism approaches 
as a basis for the discussion of our time discrete metamorphosis model and the F-convergence analysis to be 
presented in this paper. 

Flow of diffeomorphism Here, we give a very short exposition and refer to IIDGM981IBMTY051 IJMOOl 
IMTY02II for more details. Following the classical paradigm by Arnold 0Arn66llAK98l . one studies the temporal 
change of image intensities from the perspective of a family of diffeomorphisms {tp{t))te[o.i] : 5 —)■ on the 
closure of the image domain D C for d = 2,3 describing a flow, which transports image intensities along 
particle paths. In what follows, we suppose that Disa bounded domain with Lipschitz boundary. A path energy 

E[(V'(0)tG[o.i]] = / [ L[v{t),v{t)]dxdt 
Jo Jd 

is associated which each path ('i/’(f))tG[o,i] in the space of images, where v{t) = o represents 

the Eulerian velocity of the underlying flow and L is a quadratic form corresponding to a higher order elliptic 
operator. Physically, the metric 5^(t) ^(f)) = f (^)] describes the viscous dissipation in a 

multipolar fluid model as investigated by Necas and Silhavy IINv91l . From this perspective, a suitable choice 
for the viscous dissipation is given by a combination of a classical Newtonian flow and a simple multipolar 
dissipation model, namely 

L[v{t),v(t)] := ^{tie[v]f + ntr{e[v]'^) +-f\D"^v\‘^ , (1) 

wheree[u] = ^{Vv+Vv"^), m > 1 + | and A, fJ-, j > 0 (throughout this paper gradient V, divergence div, and 
higher order derivatives D™ are always evaluated with respect to the spatial variables). The first two terms of the 
integrand represent the usual dissipation density in a Newtonian fluid, whereas the third term represents a higher 
order measure for friction. Under suitable assumptions on L it is shown in IIDGM98I Theorem 2.5] that paths 
of finite energy, which connect two diffeomorphisms '0(0) = ipA and 0(1) = 0 b, are indeed one-parameter 
families of diffeomorphisms. Furthermore, for any minimizing sequence of paths a subsequence converges 
uniformly to an energy minimizing path, in particular the minimizing path solves 0(f, •) = u(f, 0(0 •)) for every 
t € [0,1], where v is the energy minimizing velocity (cf. IIDGM981 Theorem 3.1]). Given two image intensity 
functions ua,ub € L^{D), an associated geodesic path is a family of images u = {u{t) : D ®)tG[0.1] 
with u(0) = UA and m(1) = ub, which minimizes the path energy. The associated flow of images is given 
by u{t) = UA o In medical applications IIBMTY02L the diffeomorphisms represent deformations of 

anatomic reference structures described by some image ua- Thus, each diffeomorphism : D ^ for 
t G [0,1] represents a particular anatomic configuration or shape of these structures. Let us remark that this 
model is obviously invariant under rigid body motions, i.e. rigid body motions are generated by motion fields v 
with spatially constant, skew symmetric Jacobian, for which e[u] = 0 and £)"‘v = 0. 
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Metamorphosis The metamorphosis approach was first proposed by Miller and Younes OMYOlll and compre¬ 
hensively analyzed by Trouve and Younes BTYOSbl . It allows in addition for image intensity variations along 
motion paths. Conceptually and under the assumption that the family of images u is sufficiently smooth, the 
associated metric for some parameter (5 > 0 can be written as 

g{u,u) = rnin [ L[v,v] +-{u + \/u ■ v)^ dx 
J]j d 

and induces the path energy E[u] = u{t)) dt . Let = u + Vu ■ v denote the material derivative 

of u. Obviously, the same temporal change u{t) in the image intensity can be implied by different motion 
fields v{t) and different associated material derivatives i.e. u{t) = — Vu ■ v. In fact, one introduces 

a nonlinear geometric structure on the space of images by considering equivalence classes of pairs {v, ^u) as 
tangent vectors in the space of images, where such pairs are supposed to be equivalent iff they imply the same 
temporal change ii. Hence, to evaluate the metric on such tangent vectors one has to minimize over the elements 
of the equivalence class and computing a geodesic path requires to optimize both the temporal change of the 
image intensity and the motion field. Thereby, the first term L[v,v\ reflects the cost of the underlying transport 
and the term penalizes the variation of the image intensity along motion paths. 

However, typically images are not smooth and paths in image space are neither smooth in time nor in space. 
Thus, the classical notion of the material derivative u + Vu • u is not well-defined. In OTY05al Trouve and 
Younes established a suitable generalization of the above nonlinear geometric structure on L‘^{D) := K), 

which is used as the space of images, based on a proper notion of weak material derivatives. Here, we recall 
the fundamental ingr'edients of this approach. In fact, for v G L^((0,1), IL™’^(Z3, n VLq’^(D, K*^)) the 
function z G 1)) is defined as a weak material derivative of a function u G T^((0, 1),LF‘{D)) if 


n rizdxdt = — / {dtrj + div{vri))udx dt 

) Jo JD 


( 2 ) 


for rj G C“((0,1) X D). Here, IL™’^ denotes the usual Sobolev space of functions with square-integrable 
derivatives up to order m, and Wq'^ is the space of functions in with vanishing trace on the boundary. In 
terms of Riemannian manifolds, Trouve and Younes equipped the space of images LJ{D) with the following 
nonlinear structure; Let 


Nu = = (v, z) G W : J zr] + udiv{gv) dx = 0 Vp G C^{D) 

For W = n ILq’^(D,M'^)) x LJ{D) the tangent space at u G LJ{D) is defined as T^L'^iD) = 

{m} X WjNu and elements in this tangent space, which are equivalence classes, are denoted by {u, {v, z)). The 
tangent bundle is given by 

TL^{D) = U T^L^iD). 

ueL^{D) 

Furthermore, let Tr{u, {v, z)) = it be the projection onto the image manifold. Indeed, this is a weak formulation 
of the above notion of a tangent space as an equivalence class. Following the usual Riemannian manifold 
paradigm, a curve u G C'°([0,1], L^(I?)) in the space of images is called continuously differentiable, iff there is 
a continuous curve 1 1 —>■ w{t) = {v{t), z{t)) in W such that for any rj G (7“ {D) the mapping t i-G u{t)r] dx 
is continuously differentiable (denoted by u G C^([0,1], LJ{D))) and 

i z{t)T] + u{t)div{r]v{t)) dx . (3) 

In fact, for a curve t -G 7 (f) = ^it(t), (i'(f), in TLJ{D) the function z is the (weak) material derivative 

if ([^ holds for all test functions p G C“([0,1] x D) and all times t G (0,1). Furthermore, a curve u G 
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1], L‘^{D)) is defined to be regular in the space of images (denoted by u G iJ^((0,1), L^(Z3))), if there 
exists a measurable path 7 : [0,1] —?► TL‘^{D) with 7r(7) = u and bounded L^-norm in space and time, such 


that 



D 


udtrj dx dt 


n zrj + udiv{r]v) dx df 


(4) 


for all 77 G 1) X D). In fact, a continuously differentiable path u G C^([0, 1],L^{D)) is always regular, 

i.e. u G 1), L'^{D)) (cf. OTY05al Proposition 4]). Now, for a regular path u G 1), L^(£))) and 

for the quadratic form L[v, x] being coercive on H M.'^) (which can be easily verified for 

L given in ([T]) using Korn’s Lemma) one can rigorously define the path energy 



inf 




dx df. 


(5) 


In ITYOSal . Trouve and Younes proved the existence of minimizing paths for given boundary data in time. 
Adapted to our notion, they have shown that for m > 1 + | and 7 , (5 > 0 and given images ua,ub G L?{D) 
there exists a curve u G 1), L'^{D)) with u{0) = ua and u(l) = ub such that 


£\u] = inf{£’[{t] : u G 1), L^{D)), u{0) = ua, u{1) = ub} ■ 


Moreover, the infimum in 0 is attained for all t G [0,1], i.e. there exist minimizing {v, z) G 7 ( 1 ( 4 ) L\D). 

The proof relies on the observation that n Wq'‘^{D) compactly embeds into Cq'°'{D) for a < 

m- 1- f The existence of a geodesic path then follows from OTY05al Theorem 6], whereas the addendum is 
a consequence of OTY05al Theorem 2]. 


2 The variational time discretization 

In what follows, we develop a variational approach for the time discretization of geodesic paths in the metamor¬ 
phosis model. This will be based on a time discrete approximation of the above time continuous path energy 
(|^. In what follows, we suppose that 7 , <5 > 0, m > 1 + f , and define for arbitrary images u, u G L^{D) and 
for a particular energy density W a discrete energy 

>V[m, u] = min / W{D(j)) + + ^\u o cj) — u\'^ dx , ( 6 ) 

<k&AjD 0 

where A is the set of admissible deformations. Throughout this paper, we make the following assumptions with 
regard to the energy density function W : 

(Wl) W is non-negative and polyconvex, 

(W2) W{A) > /3o(det —/?! for/3o,/?i,s > 0 and every invertible matrix A with det A > 0, W{A) = 00 
for det A < 0, and 

(W3) W is sufficiently smooth and the following consistency assumptions with respect to the differential oper¬ 
ator L hold true: W{t) = 0, DW{1) = 0 and 

Furthermore, the set of admissible deformations is 

A = {4> G D) : detDf > 0 a.e. in D, </) = 1 on dD} . 
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Note that we use the symbol 1 both for the identity mapping x ^ x and the identity matrix. The first two 
assumptions ensure the existence of a minimizing deformation in (|^ and thus the well-posedness of the discrete 
energy >V[u, ft] for u, u G LF'{D). Note that BBalSlI Theorem 1] already implies the global invertibility (a.e.) 
of every (p G A because A C W^’'p{D) for a p > d. The third assumption states that the definition of W is 
consistent with the underlying dissipation described by the quadratic form L. 

Now, we consider discrete curves u = [uq, ..., uk) G (L^(D))-^^^ in image space and define a discrete 
path energy as the sum of pairwise matching functionals W evaluated on consecutive images of these discrete 
curves as follows 

K 

E^f[u] :=iT^WK_i,Ufe]. (7) 

fe=i 

We refer to IIRW14I for the introduction of such a variational time discretization on shape manifolds. Based on 
this path energy, we can define discrete geodesic paths as follows. 

Definition 2.1. Let ua,ub G L^{D) and K > 1. A discrete geodesic connecting ua and ub is a discrete curve 
in image space that minimizes 'Ek over all discrete curves u = (ug, ■ ■ ■, uk) G with ug = ua 

and Uk = ub- 


Due to the assumption (W3), the energy on the right-hand side of (|^ scales quadratically in the displacement 
(p — t, which itself is expected to scale linearly in the time step t — This already motivates the coefficient 
K in front of the discrete path energy. For the rigorous justification, we refer to the proof of Theorem 4.1 on 
the F-convergence estimates. 

In general, we want the energy density to fulfill two desirable properties: isotropy and rigid body motion 
invariance. A suitable choice for an isotropic and rigid body motion invariant energy density W in the case 
d = 2, which fulfills the assumptions (Wl-3), is given by 


= oi (tr{D(j)'^D(p)y + a 2 (det DpY + a 3 (det Dp) ® + 04 


( 8 ) 


A+/4—/4g—/xs 
+r‘s 


o ^ qrs 


with coefficients oi = 02 = ——, 1^3 - —^:;+52 

q,r > 1, which is a special case of an Ogden material. Indeed, it is possible to choose for given A, p > 0 the 
parameters q, r, s in such a way that the resulting coefficients oi, 02 and 03 are positive. Obviously, Dp 
and det Dp are invariant with respect to rotations of the observer frame. The third term of the energy density 
ensures the required response of the energy on strong compression. Both rigid body motion invariance and also 
this compression response cannot be realized with a simple quadratic energy density. For the definition of a 
corresponding energy density in the case d = 3, we refer to BCia97l Section 4.9/4.10]. 

In the discrete path energy, two opposing effects can be observed. For a given discrete curve u and (min¬ 
imizing) deformations pi,, px, the last term penalizes intensity variations along the discrete motion path 
{x, pi{x), {p20pi){x),..., {pK o ■ ■ -opYix)), whereas the first two terms penalize deviations of the (discrete) 
flow along these discrete motion paths from rigid body motions. We will see that K{uk o pk — Uk-i) reflects 
a time discrete material derivative along the above discrete motion path, whereas the first two terms represent a 
discrete dissipation density. Let us remark that minimizers of the discrete energy reversed in order are in general 
no minimizers of Ek for the reversed boundary constraint ug = ub and uk = ua- Only asymptotically in the 
limit for K -G 00 , we will obtain this symmetry based on our convergence theory below. 


3 Well-posedness of the discrete path energy and existence of discrete 
geodesics 

In this section, we will show that for images u, u G Lp{D) a minimizing deformation in the definition of 
ft] exists, which renders the definition of the discrete path energy well-posed. Furthermore, we will prove 
existence of a minimizing path u of the discrete path energy Ek and thereby establish the existence of a discrete 
geodesic. 
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Proposition 3.1 (Well-posedness of W). Under the above assumptions (Wl-2j and for u,u G Lf{D), there 
exists a deformation f G A depending on u and u such that VV[u, £t] = u, </>], where 

M, (^] := f W{D(j)) + + ^\uo (j) — dx . 

Jd ° 

Moreover, f is a diffeomorphism and 4>~^ G C^’°‘{D) for a G (0, m — 1 — 

Proof The proof proceeds in four steps. 

Step 1. Due to (Wl), we know that 0 < W := inf^g^ ft, 0] and since 1 € we have that 

1] < 00 . Consider a minimizing sequence C A with monotonously decreasing energy 

'VV^\u,u,(lf] < oo that converges to W. In particular, W = fi, < oo is an upper bound. As 

a consequence of Korn’s inequality and the Gagliardo-Nirenberg inequality for bounded domains (see 0Nir66l 
Theorem 1]), we can deduce that the minimizing sequence is bounded in W'^’^{D). Hence, due to the reflexivity 
of this space, there is a weakly convergent subsequence in again denoted by (jf, such that f 

and by the Sobolev embedding theorem we can assume uniform convergence of (f>^ f in for 

a G (0,m — 1 — |). 

Step 2. We show that the deformation f belongs to A. To this end, we will control the measure of the set 
Se = {x G D \ det D(j) < e} for sufficiently small e > 0. Indeed, by using (Wl), (W2) and Patou’s lemma, we 
obtain 


/3oe”15'.| < ^0 / {detD(l))-^dx< [ W{D<j)) dx + Pi\D\ 

Js, Js, 

<liminf/ W{D(l)^)dx +l3i\D\<W + ^i\D\ 

Js, 

and thus |S'e| < ^ which shows |S'o| = 0 and detDf > 0 a. e. on D. This implies f G A (note 

(f) G for ap > d) and due to BBalSlI Theorem 1] and f G the deformation f is injective and a 

homeomorphism. By Sard’s theorem for Holder spaces (cf. IIBHS05I 1 we additionally know that 4)~^ 

are uniformly bounded in C^’°‘{D). 

Step 3. Next, we consider the convergence of the matching functional. To this end, using the above diffeo¬ 
morphism property, we estimate 


/ \u o (j)^ — u\'^ — \u o (j) — u\'^ dx < / {\u o (j)^ — u\ + \u o (j) — u\)\u o cj)^ — u o cj)\ dx 
J D J D 

< C (IIm O (A\\l2{d) + llu O ())||i2 (£,) + ||m||l2 (£,)) \\uo(j)^ -uo fWL^^D) 

< C* (||i(||L2(_D) + ||m||l2(£))) ||m — u o 

with = (po Due to the convergence of to the identity in C^’°‘{D), we observe that the right hand 

side of the above estimate convergences to 0. To see this, we can approximate u in Lf{D) by a sequence of 
functions (uijigN and obtain 

\\U - MO < ||U - Ui\\L^(D) + \\Ui - M* Ol/''^||^2(£,) + Wu.O-tjj^ - MO'i/)-’||i2(£,) . 


The first and the third term on the right hand side converge to 0 for i —oo and fixed j, whereas the second term 
converges to zero for j —?> oo and fixed i. This establishes the convergence of the matching functional. 

Step 4. Finally, we show the lower semicontinuity for the whole functional. Let j(e) G N be such that 
[m, u,(jf] < [u, u, < W + e for all j > j (e). Furthermore, we can enlarge j (e) if necessary such 

that for all j > j{e) 


|u o (^x) — u{x)\‘^ — |{( o (j){x) — u{x)\'^ dx 


D 


< e. 
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Again using (Wl), (W2) and Fatou’s lemma, we infer 


yV^[u,u,(j)] = [ W{D(j)) + j\D^(j)\'^ + \\u o (j) — u\'^ dx 
JD 0 

< liminf f WjDd)^) + P + \ \u o — u\'^ dx + ^ < W + e + 4 , 

j^oo J d do 


which proves the claim. □ 

Next, for a given discrete path u = {uq, ... ,uk) € we define a^f^screfepaf/^ energy explicitly 

depending on a AT-tuple of deformations = {(j>i,..., (px) G as follows: 


K 




k=l 


As an immediate consequence of Proposition 3.1 there exists a vector of deformations $ G such that 
E|^[u, €>] = Eif [u]. If the images uq, ■ ■ ■ ,uk are sufficiently smooth, the corresponding system of Euler- 
Lagrange equations for (j)k is given by 

J W^A{.D(()k) ■■ DB + 2^D^(j)k : D^B + ‘^{ukO (j)^- Uk-i){Vuk o ) • 0 dx = 0 


for aWl < k < K and all test deformations B G n IEqK'^), which is a system of nonlinear 

PDFs of order 2m. Here denotes the sum over all pairwise products of two tensors. 

Before we discuss the existence of discrete geodesics, we first present the following partial result, which 
can be regarded as a counterpart of Proposition ^ . 1 [ because it establishes the existence of an energy minimizing 
vector of images u for a given vector of deformations €>. 

Proposition 3.2. Let ua,ub G and K > 2. Assume a vector $ G A^ is given. Then, there exists a 

unique u = (ug,..., Uk) G (L^(D))-^~^^ with Uq = ua, Uk = Ub such that 


Eg[u,$] 


inf E^[u,^]. 

uG(i^(r>))-^ + '-, U 0 =UA, UK = Ub 


Proof. Let u-^ = (u {,..., u^x- 1 ) *- i^))^ ^ be a minimizing sequence for the energy Ej^ [(u^, •, us), $]. 

With E|^ as a hnite upper bound for this energy along this sequence. This upper bound is obtained setting 
Uk = ^ub + {1 — ^)uA- Thanks to the estimate 

Wuih < Ihi + I O fk+l - ul\\2 + O fk+lh < (<5Ek) " + IKfc+l O </>fe + lll2 (9) 

we can deduce via induction (starting from fc = AT — 1) that u-’ is uniformly bounded in {L^{D))^~^ in¬ 
dependent of j. Thus, there exists a weakly convergent subsequence in {Lf{D))^~^ with weak limit u = 
(ui, . . .,UX-l). 

We still have to show the uniqueness of this minimizer. To this end we take into account the transformation 

rule 


/ {uk ofk- + {uk+1 o 4>k+i - UkY da^ 

JD 

= / {uk-Uk-iO(j)'i:'^)'^{detD(j)k)~^ + {uk+iO(j)k+i-Uk)'^dx 

JD 

and derive from the Euler-Lagrange equation du^ E^ [u, $] = 0 the pointwise condition 

{{uk - Uk-i O ((det D4>k)~^ o + {uk - Uk+I o '/'fe-ri)) [x) = 0 for a.e. x G D , 
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which can also be written as 


Uk+i O (j)k+iix) + {Uk-i O {x))iidet D(j)k) ^ ocl),^^{x)) 

Uk(x) = - - -^- - - ( 10 ) 

1 + (detL>(/)fe)“i o 4>k^(x) 

for a.e. x G D. This leads to a linear system of equations for {ui,..., uk-i), where evaluations at deformed 
positions are combined with evaluations at non-deformed positions, which we can consider as a block tridiago¬ 
nal operator equation. In fact, dehning for each x G D the discrete transport path 

X{x) = (Xoix), Xi(x), X 2 ix ),..., XKix)f e 


with Xo(a;) = x and Xk{x) = (j)k{Xk-i{x)) for k G {1,..., K} and the vector of associated intensity values 
C/(u,$)(a::) := {ui{Xi{x)),U 2 {X 2 ix)),... ,UK-iiXK-i{x)))'^ G (11) 


we obtain for iT > 3 and n.e. x G D n linear system of equations 


A[$](a;)(7(u, $)(x) = i?[$](a;) 


( 12 ) 


f-i 


pK-l,K-l 


is a tridiagonal matrix with 


1 


. In this case, A[€>](x) G 

( [ ]( ))fe,fe-i-i I + (det D(j)k)~^ o (j)~^(Xklx)) \ + \detD(j>k)~~^{Xk-i\x))' 

(y4[#] = + 1 , 

{detD(()k)~^{Xk-i{x)) 


(^[$](x))fc,fc_i = - 


(detD^fc) ^ o (j)y,'^{Xk{x)) 


\ + [detD(j)k) o (j)^^{Xk{x)) l + {detD(t)k) ^{Xk-i{x))' 
and i?[$](x) G is given by 


1 + (det D(j)i)~^[x) 


1 + (det D(j3K-i)~'^{XK-2{.x)) 


For any vector of regular deformations $ G A^, we recall that det D<pk > 0 for k = 1,..., Ff and $ G 
{C^(D))^. From this we deduce that for a.e. x G D the matrix A[$](x) is irreducibly diagonally dominant, 
which implies invertibility. Thus, for a\\x G D there exists a unique solution U{u, ^)(x) solving ([T^. □ 


Remark (Inherited regularity), (i) If the input images ua and mb are in L°°{D), then the images ui,..., uk-i G 
L°°{D) and they share the same upper and lower bound as the input images. This follows immediately from 
the fact that Uk{Xk{x)) can be written as a convex combination of Uk-i{Xk-i{x)) and Uk+i{Xk+i{x)) for 
fc = 1 ,..., Ff — 1 due to ( [T 0 | ). 

(ii) If the input images ua and ub are in for a < m- 1 - |, then the proof of Theorem 

shows that Uk G C^'°‘{D) for all fc = 1,..., Ff — 1. 

(iii) The intensity values along the discrete transport path X(x) depend in a unique way on the values at the two 
end points x and Xk{x) and each Uk{Xk{x)) is a weighted average of the intensities ua{x) and ub{Xb{x)), 
where the weights reflect the compression and expansion associated with the deformations along the discrete 
transport paths. 

Now, we are in the position to prove the existence of discrete geodesics making use of the existence of a 
minimizing family of deformations for the energy and a given discrete image path as a consequence of 
Proposition ITT] and the existence of an optimal discrete image path for a given family of deformations as stated 
in Proposition |3 .2 1 

Theorem 3.3 (Existence of discrete geodesics). Let ua,ub G LF'{D) and K > 2. Then there exists u G 
(L^(D))^~^ such that 

^k[{ua,u,ub)]= inf Ek[(ma,v,ub)] • 


3.2 


also 
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Proof. Let us assume that ^ with = {u{,... is a minimizing sequence of 

the discrete path energy Ek[{ua, -jUb)], where Ek is an upper bound of the discrete path energy. Due to 
Proposition |3.l[ for every there exists a family of optimal deformations , (fPj^) G with 

E^[(ua, Mb), for all . Furthermore, we can assume (by possi¬ 

bly replacing and thereby further reducing the energy) that u-’ already minimizes the discrete path energy 
E]^[(ua,v,ub),#'’] overallv e We note that due to the coercivity estimate ||Ll™(/>'il |2 < ^ 

and the Gagliardo-Nirenberg inequality the deformations are uniformly bounded in for k = 

Together with the compact embedding of into C^’°‘{D, M.^) for 0 < a < m — 1 — 

this implies that (up to the selection of another subsequence) converges to = (^i,..., fx) weakly in 

and uniformly in Following the same line of arguments as in Step 2 of the 

proof of Proposition [ 3 ^ we in addition infer that det Dfk > 0 a.e. in D for fc = 1,..., iT and thus $ S A^. 

Due to we know that the resulting images u\., which are associated with the above subsequence of defor¬ 
mations, are uniformly bounded for A: = 1,..., Ff — 1 in Lf{D)). Hence, a subsequence of converges 

weakly in L^{D) to some Uk- Finally, we deduce from the strong convergence of in that 


K 

E 


/ (ufc O (/)fc - Ufc_i)^ dx < liminf / {ui o - ui_A 

Jd ^ Jd 


‘ dx. 


Together with the weak lower semi-continuity of f i-G JjyW{D(j)) + dx we obtain with u = 

(iti,.. .,UK-i) that 

Ek[ua, u, ub] = E^[('ua, u, ub), *] 

< liminf E^[(uy4, u-t, ub), = liminf Eb-[ua,u^, wb] • 

j—foo j—foo 

This proves the claim. □ 


4 Convergence of discrete geodesic paths 

In what follows, we will study the convergence of minimizers of our discrete variational model Q for K —>■ 00 
to minimizers of the continuous model Q and thus the convergence of discrete geodesic paths to continuous 
geodesic paths. To this end, we prove F-convergence estimates for a natural extension of the discrete path 
energy. For an introduction to F-convergence, we refer to IIDal93l . 

At first, let us discuss a suitable interpolation of continuous paths. For fixed K >2 and time step size r = 
^,letffc =/cr denote the time step corresponding to a vector of images u = {uq, ... ,uk) £ For 

a vector = (^ 1 ,..., fx) £ A^ of optimal deformations resulting from the minimization in (|^, we define for 
k = 1,..., K the motion held Vk = K{(j)k — 1) and the induced transport map yk{t, x) = x + {t — tk-i)vk{x) 
with t G [tk-i,tk]- Note that yfc(ffc_i,x) = x and yk{tk,x) = (j)k{x). If one assumes that \\D(j)k — := 

supj-gB max|„|^i \{D(j){x) — l)u| < 1, then yk(t, ■) = \ + K{t — tk-i){4>k — 1) is invertible. Thus, denoting 
the inverse of yk{t, •) by Xk{t, •) one obtains the image interpolation u = UkI^, €>] with 

Uk[ii, ^]{t,x) = Uk-iixkit,x)) + K{t - tk-i){uk ofk- Uk-i){xk{t,x)) (13) 

for t G [tk-i,tk]- This interpolation represents on each interval [tk-i,tk] the blending between the images 
Uk-i = Z^ir[u, •) and Uk = ^]{tk,-) along affine transport paths 

{{t,yk{t,x)) 1 1 G [4-1,4]} 

for X G D. Based on this interpolation, a straightforward extension £k '■ T^((0,1) y. D) ^ [0,c»] of the 
discrete path energy Ek is given by 

( E^[u, $] ; if u = Uk[u, $] with u G {L^{D))^+^ and 

£k [tt] = < $ is a minimizer of E^ [u, •] over A^ 

I -foo ; else 


9 




Now, we are in the position to discuss the F-convergence estimates. The statements of the theorem are sufficient 
to prove that subsequences of discrete geodesics converge to a continuous geodesic (cf. Theorem 4.2 1 . 


Theorem 4.1 (F-convergence estimates). Under the assumptions (Wl-3}, the time discrete path energy £k F- 
converges to the time continuous path energy £ in the following sense. The estimate liminfif_>oo £k[u^] ^ 
£[u] holds for every sequence C Lf{(Q, 1) x D) with u (weakly) in T^((0,1) x D). Fur¬ 
thermore, for u G L^((0,1) X U) there exists a sequence C T^((0,1) x D) with —t u in 

L^((0,1) X D) such that the estimate limsup^_j,g^ £k\u^] < £\u] holds. 


Let us at hrst briefly outline the structure of the proof to facilitate the reading. The proof itself refers to the 
outline with corresponding paragraph headlines. To verify the liminf estimate we proceed as follows: 

(i) Reconstruction of a flow and a weak material derivative. For a sequence of images [u^, 

in L^((0,1) X D) with = (uq ,..., u^) G we consider a set of associated optimal 

matching deformations and construct the induced underlying motion held, for which the mismatch energy 
turns out to be the weak material derivative in the limit. 

(ii) Weak lower semicontinuity of the path energy. Using a priori bounds for the sequence of motion helds 
and material derivatives, we obtain weakly convergent subsequences and using a Taylor expansion of the 
energy density function W, we show a lower semicontinuity result required for the lim inf inequality. 

(iii) Identification of the limit of the material derivatives as the material derivative for the limit image se¬ 
quence. We still have to show that the pair of the weak limits of the velocity helds and the material 
derivatives is indeed an instance of a tangent vector at the limit image. The core insight is that instead of 
taking the limit in the dehning equation Q of the weak material derivative in Eulerian coordinates one 
has to use an equivalent how formulation in Lagrangian coordinates. 

(iv) Convergence of the discrete image sequences pointwise everywhere in time. In step (iii) we need that an 
image sequence with bounded path energy converges not only weakly in L^((0,1) x D), but for every 
time t G [0,1] the image sequence evaluated at that time converges already weakly in Lf{D). We use a 
trace theorem type argument to verify this. 

The proof of the lim sup estimate consists of the following steps: 

(i) Construction of the recovery sequence. The key observation is that the construction of a recovery se¬ 
quence is not based on some (time-averaged) interpolation of the given image path u G T^((0,1) x D). 
In fact, one considers for hxed K a local time averaging of an underlying motion held leading to a 
bounded path energy, and constructs from this via integration of the associated material derivative along 
the induced transport path a discrete family of images {u^, • • • , u^). 

(ii) Proof of the the lim sup inequality. The key ingredient for the proof of the lim sup inequality is the 
convexity of the total viscous dissipative functional, which we exploit based on the above construction 
of the recovery sequence via an application of Jensen’s inequality. This requires that the discrete motion 
helds are indeed dehned via local time averaging of the given continuous motion held. Furthermore, we 
again use a Taylor expansion of the energy density function W. 

(iii) Convergence of the discrete image sequences. Due to the fact that the recovery sequence of images 
{u^)kgn is dehned via integration of the material derivative and not by simple time averaging, we are 
still left to verify that converges to u in L^((0,1) x D). 
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Proof. Throughout the proof we will use a generic constant C independent of K. 

The liminf — estimate: 

(i) Reconstruction of a flow and a weak material derivative. Let C L^((0,1) x D) be any 

sequence of images that converges weakly in L^((0,1) x D) to u G L^((0,1) x D). To exclude trivial 
cases, i.e. liminf^_>,oo = oo, we may assume Sx[u^] < £ < oo for all K G N, which implies 

for = {uq,...,u^) G and an associated vector of deformations 

which is defined as a vector of (not necessarily unique) solutions of the pairwise match¬ 
ing problems (|^. Each generates on each time interval , tk) affine transport paths with motion velocity 
v^{t,y) = K{(j)^ — l){x^ {t,y)). Here, we use the notation ^ (for the sake of brevity without explicit 
reference to the sequence index K) and is the above defined pullback associated with the deformation 
on the interval [tk-iflk]- As it will be shown below in ( fTS) !, for sufficiently large K a piecewise affine recon¬ 
struction of along straight line segments from x to <f>^{x) can be performed using ( [T3] l. Thus, the difference 
quotient K {u^{4>k^x)) — u^_i{x)) is the material derivative of for all y^{t, x) with t G (tk-iflk), i-e. 


{t + s,y + SV^ if y))\^^^ = K {u^ o - u^_i) {x^{t,y)) 


is the classical material derivative of . Hence, the regularity of stated in Proposition |3. l| implies that 
fulfills the equation for the weak material derivative (|^, i.e. 


(14) 

K 


z^-ddtdx = - 


D Jo 


[ [ {dtd + div{d^d))u^ dtdx 
JD Jo 


(15) 


for all id G VLq’^((0, 1) x D) and with v^(t,y) = Vki^iV) ^ ^ [tk-iflk)- Let us remark that v^(t,-) 
vanishes on the boundary dD for t G (0,1), which corresponds to the assumption on the continuous velocity v 
in the metamorphosis model from the introduction. As a next step we show 


lim [ [ \z^'^ dtdx= lim [ \uk o fk ~'^k-i\^ ■ 

i<^^°°JdJo K^oo ^.^Jd 

Indeed, using ( [T4l l one obtains 

dtdx= j [ {{u^ {x^{t,x)))‘^ dtdx 

[ [ {{uk °(l^k-u^-i) detDy^{t,x)dtdx, 

JD Jtk — i 


(16) 




' D Jtk-i 


where Dy^{t^ x) = 1 K{t — tk-i){D(j)^(x) — 1). From the uniform bound on the energy, we deduce 


^ r 

°^k - Uk-if dx<5£ . 

k=iJD 


(17) 


Furthermore, we can estimate 

||det(l + iT(- ——1)) —< CWfk — l|lci(£>) • 

The Sobolev estimate Wf — l||ci,a(£ 3 ) CWf — for a < m — 1 — ^ and the Gagliardo-Nirenberg 

intei-polation inequality \\<j) — 1 ||w"‘. 2 (d) < 61 || 6 )™</>||j;, 2 (£)) (cf. 0Nir66l ) for f G n Wq''^{D) imply 

K 


4>k -1| 






< 


C£ 


1=1 


'LHD) - ryK ■ 


(18) 
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Together with ([T^ this proves ( [Thl l. 

(ii) Weak lower semicontinuity of the path energy. Next, from ([T^ and ( [Thl l we deduce that the material 
derivatives are uniformly bounded in T^((0,1) x D) independent of K. Thus, there exists a subsequence, 
again denoted by which converges weakly in T^((0,1) x D) to some z G ^^((0,1) x Z3) as iT —>■ oo. 

By the lower semicontinuity of the L^-norm, one achieves 



dfdx < liminf 
K—foo 



dt da;. 


Now, we will prove that there exists a velocity field v G T^((0,1), Wq (D) n such that (v, z) G 

and 

1 ^ 

f [ dxdf < liminf dx . 

Jo Jd 

The second order Taylor expansion around tk-i of the function t i-G W{t + {t — tk-i)DvjJ) att = tk gives 
W{Df^) =W(t) + ^DW{t){Dv^) + ^D^W{t){Dv^ ,Dv^) + 0{K-^\Dv^\^) 

= ^ + 0{K-^\Dv^\^) 

with {x) = K{(j)^{x) — x). The second equality follows from (W3). Then 

K 

/ W{Df^)+y\D”^f^\^dx 

k=iJD 

^^E / \{tre[v^]f+ptr{e[v^Y)+^\D^v^\^dx + cj2K jK-^\Dv^?dx. 


fc=l 


fc=l 


The last term is of order K 2 , which follows from the boundedness of the energy and by applying ( [TSl l, i.e. 

K r- K 

/ if-3|L»ufpda;<C^jnax^||(^f-l||ci(r,)E-^lk^-l|lE.2(D) <CK-^. 

fc=i k=i 


Next, for K ^ 00 the limes inferior of the remainder can be estimated as follows. We define G T^((0,1) x 
D) via v^(t, •) = vjf for t G [tk-i,tk). Due to the uniform bound of the discrete path energy is uni¬ 
formly bounded in T^((0,1), and up to the selection of a subsequence converges weakly in 

L2((0, R^) n Wq'^{D,R‘^)) to some v G ^^((0,1), K'^) n Wq'^{D,R^)) for K 00 . 

Then, by a standard weak lower semicontinuity argument we obtain 


liniinfl^ [ ^(tre[z;f]r+/rtr(e[uf]2)+7|ZZ-uf|"dx 

= liminf r [ ^{tie[v^]f + ptiie[v^]^)+y\D’^v^fdxdt 

Jq Jjy 2 

> f [ ^(tre[r;])^-f ptr(e[r;]^) + 7 dxdf . 

Jo Jd 2 

(Hi) Identification of the limit of the material derivatives as the material derivative for the limit image 
sequence. It remains to verify that we can pass to the limit in ( [T5] l for iT —> cxd with v also being the weak 
limit of in T^((0,1) x D). This will indeed imply that z is the weak material derivative for the image 
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path u and the velocity held v fulhlling Q and hence {v, z) £ Ty^L^{D). To this end, the main difficulty is to 
prove the weak continuity of [u, v) itdiv(up). In BTYOSal Theorem 2] (with the essential ingredient, which 
we actually required here, given in OTY05al Lemma 6 ]) it is shown that for the family of diffeomorphisms 
ij} : [0,1] —>■ C^{D) resulting from the transport 

Tpit,-) = (19) 

for some velocity held v £ T^(( 0 , 1 ), n Wq'^{D)) and for given initial data ' 0 ( 0 ) = 1 the integral 

formula ^ 

u{t,x) = u{0,ipt,oix)) + z{s,tjjt,s{x))ds ( 20 ) 

Jo 

for an image path u, a function z £ T^((0,1), Lp'{D)) and for a. e. x G D with ipt,s = ('0(f, O)”^) is 

equivalent to 0 . We refer to IIDGM981 Lemma 2.2] for the existence of a unique solution 0 of From ( [T^ 
we deduce that {u^, , z^) obeys 

= u^(O,0fo(a;)) + / '0f^(a;)) ds , (21) 

Jo 

where = tjj^{s, •)) with 0^ : [0,1] —?> C^{D) denoting the time discrete family of diffeomor¬ 

phisms induced by the motion held and solving 

x) = v^{t, 0^(0 x)) ( 22 ) 

for all X G D. In what follows, we will show strong convergence of 0^ to 0, for which ( [T^ holds. At hrst, 
we observe that ||2/f (0 •) ||ci,c(^) < <^(1 + l^f (0 ■)||ci.o(B)) foryf (0 a;) = x + (t - tk)v^(x) and 
t £ [tk-i,tk)- By Sard’s theorem in Holder spaces IIBHS051 and ( fTSl l we deduce that ||a:|‘^(t,')< 
(7(1 + K~^ ■)||ci.c.(£i)) for rhe inverse x^{t, ■) = Using the dehnition of vjf, the (7^’“- 

estimate for the concatenation of (7^’“-functions, and ( fTS] ) we get 

\\vkit, •)||ci.-(S) ^ hkit, •)||c.i,„(£,) (l + K-^ ||uf (0 •)||ci.o(5)) 

< C'lhf(f>')|lci.°(D) ■ 

The uniform boundedness of in L^((0,1), lU™’^(i9)) and the continuity of the embedding of 1U’”’^(Z1) 
into imply that is uniformly bounded in L^((0,1), (7^’“(Z))). Following OTY05al Lemma 7] 

(in a straightforward generalization for velocities uniformly bounded in L^((0,1), (7^’“(Z)))) one shows via 
Gronwall’s inequality that 0^ dehned in ( |2^ is uniformly bounded in L°°((0, 1),C^’°‘{D)). Finally, using 
this bound and once again the (7^’“-estimate for the concatenation of (7^’“-functions we obtain from ( |22l l the 
estimate 

||0^(0-) -0'^(s,-)||ci.c(m - 

J s 

1 

< C'(f-s)’ lh^(L')||ci.-(D) <C{t-s)^, 

which proves that 0^ is uniformly bounded in (7°’5 ([0,1], (7^’“(i9)). Thus, for some 0 with 0 < P < 
min{f, a} and up to the selection of a subsequence 0^ converges strongly in (7°’^([0,1], (7^’^(79)) to some 
0 £ (7°’5 ([0, 1 ], (7^’“(i9)) and 0 solves ( fT9| ) (cf. OTY05al Theorem 9]). The mapping (t !->■ (0^(f, ■))~^) 
which solves (|22]l backward in time, is uniformly bounded in (7°’5 ([0,1], (7^’“(i9)) (cf. OTY05al Lemma 9]). 
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Next, we obtain from ( |2T| l for functions with bounded energy the following estimate: 


< 


J z^{s,tlj^{s,x))ds^ dx 


nt-\-T 


<T\\detD{{'ip^)-^) 


K 


IL“((0,1 )xD) 
< Ct 


(sr)| 


LHD) 


ds 


(23) 


- ' Ir IIl 2 ((o,i)xd) 

for all f > 0, T > 0 with f + r < 1. The analogous estimate holds for , and replaced by u, ip, 

and z, respectively (cf. OTY05al l. From this and the uniform smoothness of ip^ and ip we deduce that for 
a subsequence (again denoted by {u^)k^'m) u^{t) u{t) weakly in L?{D) for all t G [0,1], A detailed 
verification is given in the last step of the proof below. Then, multiplying ( |2T] i with a test function ry G (£>) 
and integrating over D yields 

0=1 {t,x)ri{x)dx — f {0,ip^Q(x))ri{x) dx — f f z^{s,ip^^{x))r]{x) dxds 

JD JD ’ Jo JD 


ID 


u (t,x)ri{x)dx- / u (0,2/)p((i/'t,o) (2/))(det T'V't.o) ((V't.o) (?/))d?/ 


ID 


f [ z^is,yU{iPt)-\ymetDiPl)-\{ip{^,)-\y))dyds. 
Jo Jd 


(24) 


Based on the weak convergence of , z^ and the strong convergence of t i— [ip^(t, •)) ^ and ip^ we can 
pass to the limit in (|24ll and obtain 


0 = 


u{t,x)r]{x)dx- j u{0,y)ri{{ipt,o) (y)) (det DV't.o) ((V't.o) {y))dy 


Id 


- f [ z{s,y)ri{{iPt,s) \y)) {det DiPt,s) ^{{ipt.s) ^{y))dyds 
Jo Jd 

= / u{t,x)r]{x)dx - / u{0,ipt^o{x))r]{x)dx - / / z{s,ipt^s{x))r]{x) dxds , 

Jd Jd Jo Jd 

which shows that u and 2 fulfill ( |20l l for a. e. x G D. Since ( |20l l is equivalent to Q, this finally proves Q. 

(iv) Convergence of the discrete image sequences pointwise everywhere in time. It remains to prove that 
for a subsequence of the discrete intensity functions (again denoted by (u^)K£n) {t) u{t) weakly 

in L'^{D) for all t G [0,1]. To this end consider an arbitrary test function p G Cf°{D), t G (0,1) and r > 0 
sufficiently small (in the what follows for f = 0: f — r is replaced by t and for f = 1: f + r is replaced by t). 
Then, we obtain 


/ (t,x) — u(t,x)) ri{x) dx 

Jd 


ID 

j-t+T 


It- 


t — T J D 

.r 


Jt — T J D 

ft+r J-/ J „ 1 rt+T 


/ {u^ {t, x) — {s, x)) r]{x) — {u(t, x) — u{s, x)) ri{x) dx ds 

J D 

/ {s,x) — u{s,x)) r]{x) dxds . (25) 

Jd 


Here, f{s) ds = ^ f{s) ds is the time-averaged integral of / on (f — r, f + r). Due to the weak 


convergence of u 
AT —>■ oo. Setting 


K 


in T^((0,1) X D) the second integral on the right-hand side of (|25|) vanishes as 


ij (t,y)=v{iJ {t,y))detDip {t,y), fi{t, y) = y{ip{t,y)) det Dip{t,y) 
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we can rewrite the first term in the first integral on the right-hand side of ( |25] l and get 

rt+T 

J t — T J D 


[ - u^{s,'ip^{s,x))fi^{s,y)dyds 

J D 

[ {■u^{t,'tp^{t,y)) - u^{s,ip^{s,y))) fi^{t,y)dyds 

JD 

[ u^{s,^^{s,y)) {fi^{t,y)-fi^{s,y)) dyds. 

JD 


t — T J D 

rt+T 

± 


(26) 


t — T J D 


The second integral on the right-hand side of ( |26| ) vanishes due to the smoothness of rj and as r —> 0. 
Furthermore, using (|2^ the first integral can be estimated by 


c 


[ - u^{s,'il;^{s,y))) {t,y) dy ds 

r JD 


< sup 

sG[i—r,t+r] 

< Crs ||i7^(f,-) 


lL2(I>) 


and thus also vanishes for t —t 0. Analogous estimates apply to the remaining expression in ( |25] l replacing 
, fj^, and tjj^ by u, fj, and tp, respectively. Altogether, this proves u^{t) u{t) weakly in L^{D) for all 
fe [0,1], 

The limsup — estimate: 

(i) Construction of the recovery sequence. Consider an image curve u G L^((0,1) x D). Without any 
restriction we assume that the energy 



is bounded, where v G A^((0,1), W^'^{D) H Wq’^{D)) and z G A^((0,1) x D) are an optimal velocity field 
and a corresponding weak material derivative, respectively. Now, we define an approximate, piecewise constant 
(in time) velocity field 

for k = 1,..., K and again denoting tk = ^. 

Obviously, converges to v in Lf{{d,l),W'^’^{D)). We denote by the associated flow of dif- 
feomorphism generated by the flow equation 'tp^{t,x) = {t,tp^ {t,x)) as in ( | 2 ^ (for deduced from 
(p^ = 1 -f K~^v^) with 'ip^{0,x) = X and by ipf^ = tp^{s,{'ip^)~^{t,-)) the induced relative defor¬ 
mation from time t to time s. From this, we also obtain the underlying vector of consecutive deformations 
= {(p ^,..., p^) with tk- Following 0DGM98II we easily verify that the evolution equation 

for tp^, the uniform smoothness of tp^ and the bound on the energy 8\v\ imply that tp^ is uniformly bounded 
in C°’£[0, 1],C^'°‘{D)) (cf. the proof of the lim inf-estimate above). 

Next, the approximate discrete image path ^ • ■ •; is defined by a discrete counterpart of ( |20l l, 

namely 

uf (x) = u(0,V’£o(a;)) + / z{s,p^^^,[x))ds (27) 

Jo 

for fc = 0,..., AT. Using ( [T3] l one obtains = Uk[u^, as the requested approximation of u for given 
a: € N. 
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(ii) Proof of the the limsup inequality. At first, we verify that liinsup^_^oo £kK] < S\u]. From the 
minimizing property of Uk [u^ we deduce 

tv „ 1 

£K[u^]=-EK[n^]<KY, / WiDf^)+y\D"^f^f + -\uj^ -utifdx. 

Using the Cauchy-Schwarz inequality we derive from ( |27l i 

JD Jd 




dx 


< — 
- K 


< 


h [ [ k(s,a:)pdeti:)('0t^_^_J ^(x)dxds 

Jtk-i JD 

^ It + / l^(s>a:)Pdxds, 


where we have taken into account the estimate |1 — det g) ^(x)| < CK i, which follows from 

the uniform bound for in ([0,1],(7^’“(f5)). Furthermore, we obtain via Taylor expansion and the 
consistency assumption (W3) 

Jd 

< ^D^Wil){Dv^,Dvf) + pdx + Cj^ ±\Dv^\^dx 

= da; + l^dx. 

The definition of vjf together with Jensen’s inequality implies 

f L[v^,vj^]dx < K f [ L[v,v\dtdx . 

J ID ID J _ 1 

To estimate the remainder of the Taylor expansion we proceed as follows. At first, we obtain 

K 


\V 


K\ 


k \\C^{D) 




IT 


K\ 


1^1 


I Ww^’^iD) 


<CK f\\v{t,- 
Jo 


Iw’".2(d) df < CK 


using the Sobolev embedding theorem together with the Cauchy-Schwarz inequality and the boundedness of 
the energy E[vi\. Hence, 11^1(5) ^ CK^, which implies 

< ^i^ax^lhf|lci(c)E_^ 1^-^/ Dv{t,x)d^ dx 


< ck^ — Y 

- K 


2 K . 


K 


\Dv{t,x)\^ dtdx <CK'^ . 


k=i •' JJ J*k-1 


From these estimates we finally deduce 


£kK] < 


0 JD 


1 c 

L[v,v\ + --|zp dxdf + CK~^ + —K~^ . 
0 0 
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(Hi) Convergence of the discrete image sequences. We are still left to demonstrate that —?► m in i^((0,1) x 

D). To see this, we hrst observe that by the theorem of Arzela-Ascoli and after selection of a subsequence 
converges to '0 in 1 ], C^’“(-D)) with 0 and a < m — I — 1. From this and the quantitative control 

of the inverse of the diffeomorphisms (cf. OTY05al Lemma 9]) we deduce that its inverse, and also 
converge uniformly in x, t, and s. Thus, we get that for every t G (0,1) 

for K ^ oo. Indeed, in case of the hrst claim we argue as follows. Due to the uniform bound on 0 in L^((0,1) x 
£>) we only have to show that z(s,'ip^g(x))'^7)(s, x) dx ds converges to z(s,'ipt,s(x))'^Tt(s, x) dx ds 

for all p G C^((0, 1) X D) and q = 1,2. This is easily seen via integral transform, i.e. 


n z{s,'ftAx)yH{s,x) - z{s,'ipt,s{x)yq(^s,x)dxds 
') 


-v{s,{ipt,s) iy)){det D'ipt,s) (0t,s) {y))dyds 


where the right-hand side converges to 0 for AT —>^ 00 . The argument for u(0, •) is analogous. Hence, we can 
pass to the limit on the right-hand side of ( |27| l and achieve in analogy to the corresponding argument in the 
proof of the liminf-estimate 


(^{t,x) ^ u(O,0fo(x)) -f ^ z(s,05(^))ds^ 

(^{t,x) I-G’ u{0,tpt,o{x)) + J z{s,tpt,s{x))ds^ =u, 

where the convergence is in L^{{0, 1) x D). From this the claim follows easily. □ 

Theorem 4.2 (Convergence of discrete geodesic paths). Let ua , ub G and suppose that (Wl-3) holds. 

Furthermore, for every K G N, let be a minimizer of £k subject to (0) = ua, u^(1) = ub. Then, a 
subsequence of {u^)K&n converges weakly in L^((0,1) x D) to a minimizer of the continuous path energy S 
and the associated sequence of discrete energies converges to the minimal continuous path energy. 

Proof. The proof is standard in F-convergence theory (cf MBra02l ). Choosing = ^ub + (1 - ^)uA 
and 0 ^ = 1 we obtain an a priori bound for the discrete energy £k, which implies an a priori bound for 
z^ in L^((0, 1) X D). Using ( |2T] i, the strong convergence of 0^, and the Cauchy-Schwarz inequality we get 
that for the u^, which are associated with the minimizer of £k, the estimate < C'(IKyi||^ 2 (£,) + 

II'^^IIl2((o i)xD)) holds. From this we deduce that u^ is uniformly bounded in L°°((0,1), L^(D)). Hence, 
there exists a subsequence, again denoted by (u^)KeN, with u (weakly) in L^((0,1) x D) to some u G 

L^((0, l)x D). N ow, let us assume that there is an image path ft with g [m] < £\v\. Then, by the lim sup-estimate 
of Theorem |4. 1 [ there exists a sequence {u^)K^n with G L^((0, 1) x D) such that limsup^_^go £k[u^] < 
S\u] and together with the liminf-estimate we obtain 

£ [u] < lim inf £k[u^] < lim sup £k[u^] < £[v] , 

tT—foo 

which is a contradiction. Hence, u minimizes the continuous path energy over all admissible image paths. □ 

Remark (Inherited smoothness). Continuous solutions u of the metamorphosis model inherit for all t G (0,1) 
the regularity of the input images ua and ub (up to the Holder regularity for the exponent a). This can be 
seen as follows. For a minimizer of the continuous path energy on Lf((0, 1) x D) with u(0) G Lfi^D) and 
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u(l) S L‘^{D), Trouve and Younes give in OTYOSal Theorem 4] and BTYOSal Theorem 2] a direct representation 
of the intensity function, namely 

u{t, •) = it(0, j (det ds^ o 

for some zq G L^{D) and '0(f) = 0(f, •) the underlying flow of diffeomorphisms. Now, evaluating this equation 
for f = 1 gives 

Zo = (it(l,'0(l, •)) - u(0)) (deti:>0)“^(s) ds 

Hence, zq is as regular as ua and ub (up to the Holder regularity for the exponent a) and the same holds true 
for u{t, •) for all t G [0,1]. For discrete solutions u^, the analog statement is already given in Remark]^ 



5 Spatial discretization 

We consider a regular quadrilateral grid on a two-dimensional, rectangular image domain D consisting of cells 
Cm with m G Ic, where Iq is the index set of all cells. Based on this grid, we define the finite element space 
Vh of piecewise bilinear continuous functions (cf. BBraO?! ') and denote by the set of basis functions, 

where /at is the index set of all grid nodes Xi. Now, we investigate spatially discrete deformations '■ D ^ D 
with S (fc = 1,..., K) and spatially discrete image maps Uk '■ D —>■ K (fc = 0,..., K) with Uk G Vh 
and Uq = Ua = IhUA, Uk = Ub = IhUb- Here, Ih denotes the nodal interpolation operator. Given a 
finite element function W G Vh, we denote by 14^ = (W{xi))i^i,^ the corresponding vector of nodal values. 
Furthermore, we define a fully discrete counterpart of the so far solely time discrete path energy E^r as 
follows 


I^K,h[iUo,.-.,UK)]-.= min E’^^^[{Uo,...,UK),i<^i,...,^K)]. 

<S>k&Vl,<S>k\an=l, 

k=l,...,K 

Here, hli^o, ■ ■ ■, Uk), (4>i, ..., 4)a-)] is the discrete counterpart of and obtained by approximating the 
integrals of E^ on each cell with the Simpson quadrature rule. Here, the standard 3-point Simpson quadrature 
rule in ID is extended to 2D with 9 points using the tensor product. In our numerical experiments, this 9- 
point quadrature rule performed well. In particular, compared to lower order quadrature rules, it avoids blurring 
effects in the vicinity of image edges. Let us remark that due to the concatenation with the deformation an exact 
integration with standard quadrature rules is not possible. 

Next, we study the numerical minimization of the fully discrete energy E^ ^ for fixed (T*!,..., ^k)- For 
m G Ic the Simpson quadrature takes into account nine quadrature points. Let cc™ denote the g-th quadrature 
point in Cm and the corresponding quadrature weight for g G {0,..., 8}. Then, the entries of the weighted 
mass matrix Ma[$, tk] = jg/n basis functions being transformed via deformations $, T* 

and evaluated via quadrature are given by 


8 

Mh[<^, 4-],,, := ^ 5 ] o <l>)(x') ( 0 ^- o T-)( 4 ) . 

leic 9=0 

To evaluate the entries of this matrix numerically, we use cell-wise assembly. For m G Ic, let 0™ denote the 
basis function in the cell Cm with local index a G {0,1,2,3} and I{m, a) the global index corresponding to the 
local index a in the cell Cm, i e. Oi{m,a) = 0™ on Cm- The cell-wise assemble procedure works as follows. 
First, Ma[4), 4/] is initialized as the zero matrix. Then, for the every I G Ic and every q G {0, ..., 8 } one 
identifies the cells Cm, Cm' with $( 9 ;^) G Cm and 'l'(xg) G Cm', respectively. Finally, for all pairs of local 
indices (/3, /?') with /3, /?' G {0,1, 2, 3} one adds w^0^($(x^))0;^'(T'(x^)) to Ma[$, '^]i(m,/S),/(m',/ 9 ')- 
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Now, we are in the position to derive a linear system of equations for the vector U = (C/i,..., Uk-i) of im¬ 
ages that describes a minimizer of ^ for a hxed vector of spatially discrete deformations $ = ($i,..., 
Indeed, we can rewrite the last term in the energy ^ as follows 


E E °K) 

k^l iGic 9^0 
K 

= Y, <^k]Uk ■ Uk - t]Uk ■ Uk-I + t]Uk-i ■ Uk-i) • 

k^l 

From this, we obtain for the variation of the energy ^ with respect to the A:-th image map 

= 2 (M^[$fc, $fc] -f M^[l, 1 ]) Uk - 2 M^[$fc, l]^C7fe_i - 2 M^[$fc+i, l]i7fe+i 

for fc = 1,... ,itr — 1. In the semi-Lagrangian approach for the flow of diffeomorphisms model, a similar 
computation appears in the context of the single matching penalty with respect to the given end image (cf. 
IIBMTY05I ). For a fixed set of deformations a necessary condition for U to be a minimizer of E^ ^ is that U 
solves the block tridiagonal system of linear equations A[4>]U = R-[^], where A[$] is formed by [K — 1) x 
{K — 1) matrix blocks Ak^k' G R[$] consists of A — 1 vector blocks Rfc e with 

Afc,fc-i = — 1 ]^ , Afc^fc = ‘ffc] + 1 ], Afc_fc-i-i = —1], 

R-i = 1]^C7^ , R 2 = R 3 = ... = Rir _2 = 0 , Riv-i = , 1L]C7 b . 


The energy Tfq=o K {Wk o - C//c-iP) (x^) is convex in Uk (as a quadratic function of convex 

combinations of components of Uk) and strictly convex in Uk-i- Here, we use that the quadrature rule integrates 
affine functions exactly. Hence, E^ ^ is strictly convex in U and there is a unique minimizer U = U[$] for 
hxed This implies that A is invertible and by solving the linear system A[$]U = R[^] one computes this 
unique minimizer. Numerically, the corresponding system of linear equations (cf. line 12 of Algorithm [T]) is 
solved with a conjugate gradient method with diagonal preconditioning. 

For hxed U, the deformations $ 1 ,..., are independent of each other and thus can be updated separately. 
In the case of bilinear hnite elements, we consider only even m and replace the integrand by | At tip jn 

the quadratic form Q and correspondingly by |A™ in the energy (|^. By elliptic regularity theory, 

all of the results above directly transfer to this modihed functional. Furthermore, we use the same quadrature 
rule as before for the elastic energy and obtain the fully discrete energy 


K 


E^_J(C/o,...,C/tr),($i,...,cl>t^)] = 


E(EE w’gW{D^k{x^q)) 

k—1 l^Ic 9—0 

+7 E 


k 




iGic 9—0 


where S;i[<i), := the stiffness matrix and < 1 )^ the n-th component 

of $fc. The actual minimization of E^ ^ with respect to (the numerical solution of a simple registration 
problem) is implemented based on a step size controlled Fletcher-Reeves nonlinear conjugate gradient descent 
scheme with respect to a regularized i/t.^etric on the space of deformations OSYM07I . Thereby, the gradient 
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Data: input images Ua and Ub 

Result: approximate minimizer {Ua ~ Uq , Ui ,..., = Ub) of Ek 

1 smooth Uq = Ua and C/j = Ub with the Gaussian filter with variance cr^; 

2 for j = 1 to J do 

3 K = 2 -’; 

4 Ui^ = Ui-Uovk = Q,U...,K/2- 

5 for A: = 0 to 7^/2 — 1 do 

6 calculate $ € argmin|,gv2 l/lfe+a). ®]; 


7 

8 
9 

10 

11 

12 

13 

14 end 


t^^fe+i = C^^fe+ 2 °(l + 0.5($-l)); 

end 

repeat 

W’°^'^ = {ui,...,ui,_^)-, 

compute G argmin^g(y 2 )jr Ej^[(17A, t/s), $]; 

calculate = {U(,U^^-i) via = 24[$^]-iR[$i] ; 
until - Ui||^ < threshold; 


Algorithm 1: The alternating gradient descent scheme to compute the geodesic path. 


of the energy ^ with respect to the deformation in a direction 0 is given by 

< J(f/o,..., Uk), ..., 0 >= 

K 8 

E (E + 27E • {Mj:^Sn)^e^ 

fc=l Jg/c 9=0 n=1.2 

2 ® N 

eee < (C/fc o - Uk-i) (x^) ((Vt/fe O $,) . 0) (x^) j . 

leic 9=0 

Furthermore, we take into account a cascadic approach starting with a coarse time discretization and then 
successively refine the time discretization. In each step of this approach, we minimize the discrete path energy 
and perform a prolongation to the next finer level of the time discretization. The prolongation is based on 
the insertion of new midpoint images between every pair of consecutive images. To this end, we compute an 
optimal deformation between a pair of images and insert the middle image of the resulting warp. To improve 
the robustness of the algorithm, we additionally use a Gaussian filter with variance tr^ = | /i (color images) or 
cr^ = |/i (black-and-white images) to pre-filter the input images and damp noise, where h is the mesh size. The 
resulting alternating minimization algorithm is summarized in Algorithm[T] In the applications, it is frequently 
appropriate to ensure that deformations are not restricted too much by the Dirichlet boundary condition $ = 1 
on dD. This can practically be obtained by enlarging the computational domain and considering an extension of 
the image intensities with a constant gray or color value or by taking into account natural boundary conditions 
for the deformations. This can theoretically be justified by adding constraints on the mean deformation and the 
angular momentum. In our computations, such constraints are usually not required to avoid an unbounded rigid 
body motion component of the numerical solution. 

6 Numerical Results 

In this section, we discuss numerical results for the metamorphosis model, which are obtained with Algorithmic 
proposed in the preceding section. Besides the original model (|^, we consider a simplified model, which gives 
result of comparable visual quality with less computational effort. In the original model, we use as discussed 
above for m = 4 instead of \D'^v\'^ the term | in the quadratic form Q- The simplified model is associated 
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K = i, 

using (^, 
with 8 — 10“^, 

A= i, 

<i = r=l 
s= i,7 = 10-5 





fA 

H 


K = A, 
using {29] with 
7 = 10-5, 

5 = 10-1 





-^\\j 



K = 16, 
using {29] with 
7 = 10-5, 

S = 10-1 


p7..A I p7,.A. I p-t-A p->A 



f .x f , A f f f .y\. 

v_y \_y 





Figure 1; Metamorphosis for two slices of a MRT data set of a human brain (data courtesy of H. Urbach, 
Neuroradiology, University Hospital Bonn). We compare the original model (first row) with the simplified 
model and K = A (second row), K = IQ (third to fifth row). 


with the quadratic form 

L[v{t),v{t)] := Dv : Dv + jAv ■ Av , 

where 7 > 0. A choice for the discrete energy, which is consistent with this quadratic form, is given by 

yy [u, m] = min / D(j) : Dcj) + -fAcj) ■ Acj) + \ \u o (j) — u\'^ dx . 

^ Jd ^ 


(28) 


(29) 


In fact, this retrieves a very basic model for the registration of the two images u and u consisting of a simple 
thin plate spline regularization and the most basic fidelity term (cf. OMF03I 1. 

Let us emphasize that in the spatially continuous setting both the existence theory and the F-convergence 
result require the full set of assumptions. In particular, the definitions ( |28| ) for L[-, •] and ( |29| ) for W (contrary to 
the full model with W proposed in (|^) do not comply with (W2) and (W3). However, in case of the definition 
( [29l l, the regularization term of the deformation energy W is quadratic and enables a significant speedup of the 
algorithm compared to the theoretically justified fully nonlinear model. We compare both models in our first 
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K = A, 
using {29} with 
7 = 10-3, 

S = 10-3 


K = 16, 
using {29} with 
7 = 10-3, 

5 = 10-3 


999999 

<W <y <y <y 

Figure 2: Metamorphosis between two faces from female portrait paintings. 


example and use the simplified model in all other applications. The parameter threshold is set to 10“® in the 
algorithm. 

Figure [^depicts a discrete geodesic path obtained with the full model (with parameters K = A, S = 10“^, 
A = 1, fj, = q = r = f,s = I and 7 = 10“®) and with the simplified model (with parameters K G {4,16}, 
7 = 10“®, 6 = 10“^), where ua and ub are different slices of a 3D magnetic resonance tomography of a 
human brain. 

Figure 1^ shows a geodesic path between two faces from female portrait painting^ computed with the sim¬ 
plified model with parameters 7 = 10“® and <5 = 10“^. The local contributions 'E^[{Uk-i,Uk),^k] for 
k = 1,... ,K of the total energy and its components are shown in FigureNote that the method seems to 
prefer an approximate equidistribution of the total path energy in time. 

Finally, we consider time discrete geodesic paths in the space of color images. To this end, we take into 
account a straightforward generalization of the model for scalar (gray) valued image maps to vector-valued 
image maps. One can even enhance the model with further channels. Such additional channels can represent 
segmented regions of the images, which one would like to ensure to be properly matched by transport and 
not by blending of intensities. The only required modification of the method is that \uk+i o (j)k+i — Uk\ is 
now the Euclidean norm of the (extended) color vector. As an application, we considered the metamorphosis 
between two self-portraits by van Gogh (see Figure[^[^ Since the background colors of both self-portraits differ 
considerably in the RGB color space, we adjusted the background color of one of the images (i.e. replacing 

'first painting by A. Kauffmann (public domain, see http://commons.wikimedia. 0 rg/wiki/File:Angelika_ 
Kauffmann_-_Self_Portrait_~_1784 . jpgl, second painting by R. Peale (GFDL, see http://en.wikipedia.org/ 
wiki/File:Mary_Denison.jpgl 

^both paintings by V. van (jogh (public domain, http : // en. wikipedia .org/wiki/File: SelbstPortrait_VG2 . jpg 
http://upload.wikimedia.org/wikipedia/commons/7/71/Vincent_Willem_van_Gogh_102.jpgl 
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iiiiliiili 

Figure 3: Energy contributions of the regularization functional ■ D^k + 7A<i>fe • dec (red) and the 

matching functional ^ \Uk o — Uk-i\'^ dec (green) for the discrete geodesic path in Figurej^with K = A 
(left) and K = 16 (right). 


0.00012 




Figure 5; Original van Gogh self-portraits and the background modulated input image 

together with the associated fourth channel segmentations u\ and u^. 


Ms by in Figure]^. In this application, a fourth (segmentation) channel is used to ensure the proper 

matching of the ears and the clothing. The time-discrete geodesic path for 
the van Gogh self-portraits is shown in Figure]^ for K = 8 along with the 
temporal change of the fourth channel. Again, we used the simplified model 
with parameters 7 = 10“^ and 6 = 10“^. Figuredepicts the pullback 


.RGB . 


o < 1 ) along the flow induced deformation $ = ° ^k-i o ... o <i)i 


corresponding to the geodesics in Figure]^ Finally, Figure |7] visualizes the 
deformations and the corresponding accumulated weak material derivative 
along the discrete geodesic path. The color wheel on the lower left in the 
first row indicates both the direction and the magnitude of the discrete ve¬ 
locities — !)■ Obviously, the motion field is not constant in time. 

Furthermore, to visualize the change of the image intensity along motion 
paths, the accumulated weak material derivative Zi {I = 1 ,... , 8 ) with 
Zi = O - Uk-i) O Xk -1 using the notation 0 is plotted using an equal rescaling for all 

1 . 



RGB o $ of 
image map along the path. 


Figure 4: Pullback it| 


7 Conclusions and outlook 

We have developed a robust and effective time discrete approximation for the metamorphosis approach to com¬ 
pute shortest paths in the space of images. Thereby, the underlying discrete path energy is a sum of classical 
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Figure 6; Metamorphosis between two “van Gogh self-portraits” using the energy ( |29l l for X = 8 and S = 10 ^ 
including the fourth (segmentation) channel (bottom row). 
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Figure 7: Discrete motion fields K{^k — 1) (first row) and accumulated weak material derivative Zi (second 
row) for fc = 1,..., 9. 


image matching functionals. The approach allows for edge type singularities in the input images. We have 
proven existence of minimizers of the discrete path energy and convergence of minimizing discrete paths to a 
continuous path, which minimizes the continuous path energy. This analysis is based on a combination of the 
variational perspective of (discrete) geodesics as minimizers of the continuous © and discrete path energy 0 , 
respectively, with the continuous (([T^, ( |20| l) and discrete flow perspective ((|22ll, (|2T|l). In particular, this com¬ 
bination is the basis of a compensated compactness argument for the weak material derivative. Indeed, using 
the flow perspective ( |20l l, we are able to compensate for the loss of compactness in time, when trying to pass 
to the limit in the weak definition of the discrete material derivative Using a finite element ansatz for the 
spatial discretization, a numerical algorithm has been presented to compute discrete geodesic paths. Qualitative 
properties of the algorithm are discussed for three different examples including an application to multi-channel 
images. Particularly interesting future research directions are 

- the use of duality techniques in PDF constraint optimization to derive a Newton type scheme for the 
simultaneous optimization of the set of deformations and the set of images associated with the discrete 
path, 

- a full-fledged discrete geodesic calculus based on the general procedure developed in 1IRW131|RW14I and 
including a discrete logarithmic map, a discrete exponential map, and a discrete parallel transport, and 

- a concept for discrete geodesic regression and geometric, statistical analysis in the space of images. 

Furthermore, the close connection to optimal transportation offers interesting perspectives, which should be 
exploited. 
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